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ABSTRACT 

We present a detailed numerical study of the Gough & Mclntyre model for the so- 
lar tachocline. This model explains the uniformity of the rotation profile observed in 
the bulk of the radiative zone by the presence of a large-scale primordial magnetic 
field confined below the tachocline by flows originating from within the convection 
zone. We attribute the failure of previous numerical attempts at reproducing even 
qualitatively Gough & Mclntyre's idea to the use of boundary conditions which in- 
appropriately model the radiative-convective interface. We emphasize the key role of 
flows downwclling from the convection zone in confining the assumed internal field. We 
carefully select the range of parameters used in the simulations to guarantee a faithful 
representation of the hierarchy of expected lengthscales. We then present, for the first 
time, a fully nonlinear and self-consistent numerical solution of the Gough & Mclntyre 
model which qualitatively satisfies the following set of observational constraints: (i) 
the quenching of the large-scale differential rotation below the tachocline - including 
in the polar regions - as seen by helioseismology (ii) the confinement of the large-scale 
meridional flows to the uppermost layers of the radiative zone as required by observed 
light element abundances and suggested by helioseismic sound-speed data. 
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1 INTRODUCTION 

The presence of the tachocline, a thin shear layer located at 
the interface between the radiative and convective regions 
of the Sun, was established two decades ago (Christensen- 
Dalsgaard & Schou, 1988; Kosovichev, 1988; Brown et al 
1989; Dziembowski et al. 1989) but its modus operandi still 
remains mysterious. 

Anisotropic turbulent stresses associated with rotation- 
ally constrained eddies are thought to maintain the differ- 
ential rotation profile observed within the convection zone: 



a. 



-) ~ O eq (r)(l — 0,2(1") cos 2 6 — 04,(1") cos 4 6) , (1) 



where for example at r = 0.75-Rq Q, cq /2n = 463nHz, 
a 2 = 0.17 and a 4 = 0.08 (Schou et al. 1998; Gough 2007). 
However, as shown by Spiegel & Zahn (1992) (SZ92 here- 
after), the mere reduction in the amplitude of these stresses 
naturally expected to occur across the radiative-convective 
interface (at r cz = 0.713i?o) cannot explain the transition to 
near-uniform rotation below. It was later argued by Gough 
& Mclntyre (1998) (GM98 hereafter) that only long-range 
stresses can explain the suppression of the rotational shear 
in the bulk of the radiative zone and in addition main- 
tain the angular velocity of the interior (observed to be 
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n rz /27r ~ 430nHz) close to that of the surface despite the 
global spin-down induced by magnetic-braking. 

Two competing theories for these presumed long- 
range stresses have been investigated: purely hydrodynamic 
stresses in the form of gravity waves (see the review by Zahn, 
2007) and hydromagnetic stresses (see the review by Ga- 
raud, 2007). This paper focuses on the latter mechanism 
only, although the dynamical balance in the Sun could ar- 
guably involve a combination of the two. 

It has long been known that even a very small primor- 
dial magnetic field embedded in the radiative zone could 
in principle impose uniform rotation (Ferraro, 1937; Mestel, 
1953; Mestel & Weiss, 1987). Ferraro's isorotation law, 



B ■ Vfl = 



(2) 



valid in the limit of negligible dissipation and for steady- 
state, axisymmetric flows, is usually stated as the angular 
velocity must be constant on magnetic field lines. Thus in its 
simplest form, Ferraro's law predicts the possibility of uni- 
form rotation for the radiative zone provided the magnetic 
field is entirely confined beneath the radiative-convective 
interface (Riidiger & Kitchatinov 1997). On the other hand, 
any field line directly connected with the convection zone 
promotes the propagation of the rotational shear into the 
radiative zone (MacGregor & Charbonneau 1999), inducing 
what will be referred to from here on as a "differentially ro- 
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tating Ferraro state" . Hence, field confinement is the key to 
the existence of a tachocline. 

GM98 were the first to address the question of how 
the presumed primordial field could indeed be confined 
within the radiative zone, and proposed that meridional 
flows driven by Coriolis forces in the convection zone and 
burrowing downward would interact nonlinearly with the 
underlying magnetic field lines, bending them towards the 
horizontal in the tachocline region, thus effectively suppress- 
ing direct radial Alfvenic transport between the convection 
zone and the radiative zone. Conveniently, the same merid- 
ional flows can also be held responsible for mixing light el- 
ements such as Li and Be between the convection zone and 
their respective nuclear-burning regions, reducing He set- 
tling (Elliott & Gough, 1999) as well as providing weak 
but sufficient angular-momentum transport to adjust the 
mean rotation rates of the convective and radiative zones 
continuously throughout the spin-down phase. The seminal 
boundary-layer analysis of the dynamics of the tachocline 
presented by GM98 validated the plausibility of this theory, 
although many of their simplifying scaling assumptions re- 
main to be confirmed through the direct numerical solution 
of the governing equations. 

This paper first briefly reviews existing work and then 
presents new results on the laminar tachocline dynamics ac- 
cording to GM98. We begin by discussing past attempts 
at implementing their model numerically in Section [5] In 
particular, we argue that the failure of previous numerical 
studies in reproducing field confinement can be explained by 
the selection of inappropriate boundary conditions. A new 
numerical algorithm is then presented in Section [3] and used 
in Section |4]to revisit the idea proposed by GM98 with much 
more success. We discuss future prospects in Section [5] 



2 DISCUSSION OF PREVIOUS NUMERICAL 
MODELS IN THE LIGHT OF GOUGH & 
MCINTYRE'S ORIGINAL IDEA 

GM98 argue that the radiative interior should be divided 
into three dynamically distinct zones including (from the 
base of the convection zone downward) (1) a more-or- 
less magnetic-free region in thermal-wind balance, well- 
ventilated by meridional flows originating from the con- 
vection zone, which can be thought of as the bulk of the 
tachocline, (2) a very thin magnetic advection-diffusion layer 
where the tachocline flows and the underlying field interact 
nonlinearly to confine one another and (3) a magnetically 
dominated, near-uniformly rotating interior. 

The nonlinear nature and geometric complexity of the 
problem precludes analytical solutions, and all existing stud- 
ies since the original work of GM98 have been numerical (see 
the review by Garaud, 2007). Among these, only two include 
all the nonlinear terms required (i.e. the advection terms 
in the momentum and the magnetic induction equations, 
and the Lorentz force in the momentum equation) to repre- 
sent the nonlinear magnetohydrodynamics of the tachocline 
"correctly": Garaud (2002) - hereafter Paper I - and Brun 
& Zahn (2006) - hereafter BZ06. Surprisingly, neither have 
been able to find evidence for the kind of dynamical balance 
proposed by GM98, and the time is therefore ripe to take a 
step back and discuss why. 



2.1 "Failure" of previous numerical models 

Paper I and BZ06 are two studies which model magneto- 
hydrodynamic perturbations induced in the solar radiative 
zone by the differentially rotating convection zone and by an 
assumed primordial magnetic field. Both papers focus on the 
issues of field confinement and the suppression of differential 
rotation below the tachocline. 

Paper I presents steady-state, axially symmetric solu- 
tions of an incompressible and unstratified version of the 
GM98 model. BZ06 on the other hand use a time-dependent, 
three-dimensional algorithm based on the ASH code (Glatz- 
maier 1984; Clune et al 1999; Miesch et al. 2000; Brun et al. 
2004), and solve the complete set of anelastic MHD equa- 
tions. 

Both studies otherwise consider a similar computational 
domain, namely a spherical shell which spans the region be- 
tween the base of the convection zone (at r = r out — r cz ) and 
an inner sphere (at r — r- m ). The boundary conditions are 
also essentially similar: the "outer" boundary, which effec- 
tively models the radiative-convective interface, is in both 
cases assumed to be impermeable and rotating differentially 
with an angular velocity profile given by equation JTJ; the 
magnetic field is matched onto a potential field. 

In Paper I, the numerical solutions only show confine- 
ment of the magnetic field lines in the equatorial regions 
for low enough diffusivities, occasionally at mid-latitudes 
for lower field strength, but never in the polar regions (see 
Figure 11 of Paper I for example). It is commonly argued 
that the failure of Paper I to find fully-confined magnetic 
field solutions stems entirely from the simplified nature of 
the model equations (incompressible, unstratified): indeed, 
within these assumptions the meridional flows are gener- 
ated only by Ekman-Hartmann pumping on the boundaries, 
and their amplitude scales with the diffusivities in a way 
which always maintains the magnetic Reynolds number be- 
low unity. As a result, these flows are unable to confine the 
magnetic field and the differential rotation imposed at the 
upper boundary of the domain persists within a large part 
of the radiative zone. 

The three-dimensional, time-dependent solutions of 
BZ06 naturally depend on the initial magnetic field con- 
figuration selected. Various cases with an initial field more- 
or-less deeply confined are discussed. Against expectations, 
BZ06 find that regardless of the initial conditions, the field 
lines always spread out and eventually overlap with the con- 
vection zone, permitting the propagation of the differential 
rotation into the radiative interior. Interestingly, the tran- 
sient state of the system - prior to any field line connect- 
ing with the convection zone - qualitatively looks in many 
ways similar to the well-known GM98 picture, although it 
is clearly not a steady state. The "tachocline" thus formed 
slowly narrows with time until field lines spread into the con- 
vection zone, at which point a differentially rotating Ferraro 
state is rapidly established. One is forced to conclude in one 
of three ways: (1) for unspecified primordial reasons, the ini- 
tial magnetic field was more deeply embedded in the radia- 
tive zone and the currently observed tachocline is merely a 
transient phase; (2) the tachocline is indeed in a steady-state 
and some of the assumptions made in the way have been 
wrong; (3) the parameter regime studied by BZ06 (where 
diffusion coefficients are artificially increased by many or- 
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ders of magnitude to satisfy numerical constraints) does not 
appropriately reproduce the solar interior dynamics. 

One should be uncomfortable in selecting the first op- 
tion, as it would place very strong and unlikely constraints 
on the initial field conditions to provide just the right struc- 
ture for today's tachocline. One could naively tend to favour 
the third option, but BZ06 a priori took care to select a 
range of parameters for which all expected boundary layer 
thicknesses are small enough and in the same hierarchical 
order as those of the model proposed by GM98. The plot 
thickens while we are left with the uneasy task of reconsid- 
ering the key features of either the numerical models or of 
the GM98 model (both, perhaps). 

2.2 The source of the problem 

At this point it is worth discussing one of the more deli- 
cate aspects of the GM98 model, namely the exact mech- 
anism by which these pivotal meridional flows are thought 
to be generated. This point has been a source of confusion 
and debate, but is clearly crucial to a better understanding 
of the tachocline dynamics. In the original work of GM98, 
the principal clue to the nature of the flows can be found 
in the sentence: "Turbulent stresses in the convection zone 
induce (through Coriolis effects) a meridional circulation, 
causing the gas from the convection zone to burrow down- 
wards ..." . The flows considered by GM98 do not originate 
from within the tachocline and can therefore not be ap- 
propriately modelled by a numerical scheme in which the 
radiative-convective interface is assumed to be impermeable 
(as in Paper I and BZ06). While this conclusion seems ob- 
vious in hindsight, the physics of the problem are actually 
rather subtle and deserve clarification. In what follows, we 
discuss the issue in more detail and present a unified view 
of the results of previous works on the subject. 

Spiegel & Zahn were the first to study the dynam- 
ics of the newly-discovered tachocline (SZ92). They con- 
sidered a non-magnetic radiative zone only, and imposed 
a latitudinally-varying rotation profile at the radiative- 
convective interface. They performed two distinct calcu- 
lations. The first looked at the time-dependent evolution 
of the angular velocity profile in the radiative zone un- 
der such forcing, assuming isotropic viscous stresses. The 
second sought steady-state "tachocline" solutions assuming 
anisotropic turbulent stresses. Since the latter did not specif- 
ically address the question of the meridional flows, we focus 
here on the results of the time-dependent calculation, which 
can be interpreted in the following way. The differential ro- 
tation imposed by the convection zone to the top of the 
radiative zone inevitably induces some degree of shear along 
the rotation axis if the radiative zone is originally in a state 
of uniform rotation. This generates a thermal gradient in 
the latitudinal direction by way of the thermal-wind equa- 
tion (see GM98, equation (1) for example). Meridional flows 
are then required to balance this thermal gradient when the 
system is in thermal equilibrium. These flows burrow into 
the radiative zone, advecting angular momentum thereby 
helping the propagation of the shear further down. The re- 
sults of SZ92 imply that the system continues evolving in 
this fashion until the radiative zone has achieved complete 
thermal and dynamical equilibrium. 

The characteristic amplitude of the time-dependent 



flows associated with this thermo-dynamical relaxation pro- 
cess is, by way of the assumptions listed above, only depen- 
dent on the (evolving) local differential rotation, stratifica- 
tion and thermal conductivity. Their turnover timescale is 
calculated to be of the order of the local Eddington- Sweet 
timescale, which is short in the initial relaxation stages and 
steadily increases as the system evolves. Crucially, this re- 
sult was derived independently of any boundary condition 
on the meridional flows. It is therefore correct to think of the 
induced transient flows as being driven by baroclinic stresses 
from within the tachocline rather than downwelling from the 
convection zone. In fact, they exist even if the radiative- 
convective interface is assumed to be impermeable. This fact 
is probably at the origin of the impermeable boundary con- 
ditions selected by Paper I and by BZ06, in spite the obvious 
contradiction with GM98's intent. 

However, it is vital to remember that this transient 
phase and its associated flows both end once the system 
achieves thermal and dynamical equilibrium. The only flows 
remaining in the final steady-state are driven, within tiny 
boundary layers, by the thermo(-magneto)-viscous stresses 
required to match the bulk equilibrium solutions to the ap- 
plied boundary conditions. Thus, as expected from any ellip- 
tic problem, the nature of the boundary conditions selected 
entirely controls the steady-state solutions. This is clearly il- 
lustrated by the work of Garaud & Brummell (2008) (GB08 
hereafter), which complements that of SZ92 by calculating 
the spatial properties as well as characteristic amplitudes 
of steady-state meridional flows in the radiative zone, as in- 
duced by various kinds of forcing applied at the radiative- 
convective interface. 

GB08 showed that in the non-magnetic case, the in- 
duced steady-state flows can always be viewed - at least in 
the linear sense, and for solar parameters - as the sum of 
two "modes" : a viscously-dominated "Ekman mode" , which 
very rapidly decays away from the interface (on an Ekman 
lengthscale) , and a "thermo- viscous mode" which can essen- 
tially span the entire radiative zone when not hindered by 
a magnetic field. Since the Ekman mode decays too rapidly 
to have any effect on the tachocline dynamics, the flows dis- 
cussed by GM98, which are thought to ventilate the bulk of 
the tachocline, should be identified with the thermo-viscous 
mode. This view is perfectly consistent with GM98's anal- 
ysis, since the thermo-viscous mode is in thermal-wind bal- 
ance (GB08). 

The respective amplitudes with which the Ekman and 
thermo-viscous modes are driven depend on the nature, am- 
plitude and spatial structure of the forcing. GB08 showed 
that when the interface is impermeable, as assumed in Pa- 
per I and by BZ06, then the amplitude of the thermo-viscous 
mode is negligible for microscopic solar values of the diffu- 
sivities. As a result, the magnetic Reynolds number of these 
tachocline flows is much smaller than unity, which straight- 
forwardly explains why field confinement eluded these two 
previous studies. 

When flows are pumped directly through the interface 
by stresses within the convection zone (i.e. the radial compo- 
nent of the flows is non-zero at the radiative-convective in- 
terface), GB08 find that the amplitude of the thermo-viscous 
mode can be much higher, in which case the tachocline flows 
may be expected to confine the field in a scenario qualita- 
tively similar to the one proposed by GM98. The quantita- 
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tive detailed analysis must be done numerically; this is the 
purpose of the present study. 



3 THE NUMERICAL MODEL 

We now revisit the idea proposed by GM98 in the light of 
the previous discussion. A new numerical algorithm has been 
constructed specifically for this purpose. We first present, for 
the sake of completeness, a derivation of our model equations 
in Section T3. II Since our computational domain is limited to 
a spherical shell within the radiative zone, we then present 
and discuss in detail all the boundary conditions applied in 
Section [3.21 The relevant non-dimensional parameters, and 
expected boundary layers of the problem are discussed in 
Section 13.31 The numerical method, as well as other numer- 
ical constraints, are finally outlined in Section 13.41 



3.1 The governing equations 

The Standard Solar Model (SSM hereafter) provides an ac- 
curate representation of the thermal structure and composi- 
tion of the interior of a hypothetical spherically symmetric, 
non-rotating, non-magnetic Sun. The excellent match be- 
tween the SSM and helioseismic observations confirms that 
the likely internal dynamics of the Sun only induce weak 
deviations in the background state thermodynamical quan- 
tities. The large-scale magnetohydrodynamics of the solar 
interior, thought to be responsible for the maintenance of the 
peculiar observed rotation profile, can therefore be thought 
of as perturbations upon the SSM. 

The quasi-steady SSM equations (which are valid on 
timescales much shorter than the nuclear burning timescale) 
reduce to the following set of four equations within the solar 
radiative zone (but outside of the nuclear burning core): 

-Vp = pV$ , 
V-(fcVT) = , 
P = P(P,T) , 

V 2 ¥ = 4irGp . (3) 

These equations describe hydrostatic equilibrium, thermal 
equilibrium, the equation of state and finally the Poisson 
equation for the gravitational potential respectively. Here, 
p, T, and p are the standard thermodynamical variables, k is 
the thermal conductivity (which typically depends on T and 
p) , G is the gravitational constant and $ is the gravitational 
potential. The equation of state itself is well-approximated 
by that of a perfect gas in this region of the Sun. The solution 
to the set of equations @ for the present-day Sun can be 
inferred from model S of Christensen-Dalsgaard et al. (1991) 
for example (see Appendix A). 

We then consider perturbations on this background 
equilibrium caused by meridional and azimuthal flows, as 
well as magnetic fields. For simplicity, we restrict our study 
to axisymmetric systems. Moreover, since the observed inter- 
nal rotation profile of the Sun has remained approximately 
constant over the past decade, we boldly postulate that the 
Sun is in fact in a state of quasi-steady dynamical and ther- 
mal balance. The general set of equations governing the ax- 



isymmetric, quasi-steady perturbations are then 

pu-Vu = - Vp - pV$ + j x B + VTi , 

V-(pu) =0 , 

pTu ■ Vs = V-(fcVT) , 

P = P(P,T) , 

V 2 <!> = 4-xGp , 

Vx(u x B) = Vx(jiVx.B) , 

V B = , (4) 

where u represents the axisymmetric flow velocity, B the 
magnetic field and j — X7xB/4n is the electric current, II 
is the viscous stress tensor, s is the specific entropy and rj is 
the magnetic diffusivity. These equations respectively char- 
acterise the conservation of momentum, mass, and thermal 
energy, the equation of state and the Poisson equation, con- 
servation of magnetic flux and finally, the solenoidal condi- 
tion. 

Thermodynamical perturbations are from here on de- 
noted with tildes, as for example in p = p + p. The sys- 
tem of equations (j4| can be linearised in the thermodynam- 
ical quantities around the background state equations ((3} to 
yield 

pu-X7u = -Vp - pV$ - p"V$ + j x B + vn , 

V-(pu) = , 

pTu- Vs = V-(feVT) , 

l = l + t 
P p T ' 

V 2 <§ = 4-xGp , 

Vx(it x B) = Vx^VxB) , 

V-B = 0, (5) 

where we approximated the equation of state with that of 
a perfect gas, and neglected perturbations to the chemical 
composition. Note that the latter approximation was chosen 
for the sake of simplicity: Wood & Mclntyre (2007) showed 
that the effects of compositional gradients (notably Helium) 
could play an important role in the tachocline dynamics. The 
background magnetic diffusivity rj, kinematic viscosity V and 
thermal conductivity k are calculated from model S using 
the expressions provided by Gough (2007) (see Appendix 
A). 

Note that the nonlinearities in the quantities relating 
to the magnetohydrodynamics of the interior (flow velocities 
and magnetic field) have not yet been tampered witlQ. The 
next step is to use a spherical coordinate system (r, 9, (j>) 
and move to a rotating frame of reference, by considering 
that u = u + u, with u — r sin O^le^ and u = (u r , ue,u$)- 
The value of Q adopted is uniquely defined by requiring that 
the total angular momentum of the convection zone be null 
in the rotating frame (Oilman et al, 1989): if tt C z(6,r) is 
assumed to be independent of radius, and given by equation 

Also note that in principle, the thermal energy equation should 
contain an additional term, namely fcVT. This term is included 
in the numerical algorithm for the sake of completeness, but in 
practise has little influence on the result (at least in the case of 
the solar radiative zone). It will not be discussed in this paper for 
the sake of simplicity. 
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CE}, then 

n -°-('-T-£)- (6) 

The equation for the conservation of momentum then be- 
comes 

pu ■ Vm + 2pfl x u + pfl x CI x r 

= -Vp-pV$-pV$ + j x B + v-n , (7) 



while the assumption of axial symmetry implies that the 
equations for the conservations of mass, thermal energy and 
magnetic flux are merely changed by replacing u with ii. 

It is well known that the baroclinic deformation of 
the background state caused by the centrifugal force drives 
meridional flows, albeit very slow ones since the Sun is a slow 
rotator (Eddington, 1925; Sweet, 1950). Since the global 
turnover timescale of this Eddington-Sweet circulation is 
calculated to be orders of magnitude longer than the age of 
the Sun, it cannot play any role in the angular momentum 
or chemical mixing processes. In a steady-state calculation, 
however, these flows play an artificially important role and 
it is important to suppress them. This can be done by us- 
ing the following momentum conservation equation instead 
(SZ92): 

pu ■ Vu + 2pTi x u = -Vp - pV$ + j x B + V-n , (8) 

where Cowling's approximation was also used to justify ne- 
glecting the perturbation of the gravitational potential re- 
lated to density and temperature perturbations caused by 
the remaining Coriolis and Lorentz forces. 

Finally, we assume that the meridional flow velocities 
remain small (SZ92; GM98): we neglect quadratic terms in 
u r and He, but retain nonlinearities in u$ to allow for sig- 
nificant differential rotation. 

To summarise, the model equations we consider are 

pu ■ Vw + ZpTi x u = -Vp - pV$ + j x B + v-n , 

V-(pu) =0 , 

pTii- Vs = V-(AVT) , 

P = E + t 
P p T 

Vx(u x B) = Vx(r/VxB) , 

V-B = , (9) 

where the quantities solved for are the three components 
of u (as defined above), the three components of B = 
(B r , Be, B^), and the three thermodynamical variables p, 
p and T, and where quadratic terms in u 2 ., Hg and u r ue 
are neglected (see Appendix B for the full set of equations 
expanded in a spherical coordinate system). 

This set of equations is equivalent to the one used by 
GM98 before they apply their boundary-layer analysis. It 
also reduces to the one used by SZ92 in the non-magnetic 
case, provided the following further approximations are ap- 
plied: the Boussinesq approximation, assuming that Q, <C SI 
and that ue u r . Finally, it is a steady-state and axisym- 
metric version of the equations used by BZ06. 

For comparison, note that Riidiger & Kitchatinov 
(1997), MacGregor & Charbonneau (1999), and Kitchati- 
nov & Riidiger (2006) neglect meridional flows entirely in 




r in-0-35Rg, r out -0.7RQ Rq 



Figure 1. The computational domain is limited to the spherical 
shell shown here in grey. 

their calculations of the structure of the radiative interioi0 
Sule, Arlt & Riidiger (2005) include meridional flows but 
assume that the poloidal field is fixed, so that the interac- 
tion between the field and the flows through the Lorentz 
force in their paper is fundamentally linear (and therefore 
inappropriate) . 

3.2 Domain boundaries and boundary conditions 

3.2.1 Computational domain 

The geometry of the computational domain is shown in Fig- 
ure [T] The upper boundary is selected to be slightly below 
the base of the convection zone, at r out = 0.7Rq. The lower 
boundary is at ri n = 0.35-Rq as in BZ06. The inner core 
(for r < ri n ) can be viewed as a solid metallic sphere, within 
which a permanent magnetic dipole is maintained. The selec- 
tion of the position of the lower boundary was found to have 
little influence on the tachocline dynamics for low enough 
values of the diffusivities. 

3.2.2 Model boundary conditions 

The selection of adequate boundary conditions for this 
model is a very delicate task. Given the elliptic and non- 
linear nature of the system considered, the existence of so- 
lutions is not guaranteed, and solutions (when they exist) 
are entirely controlled by the boundary conditions applied. 

Our aim is to propose a set of boundary conditions at 
the upper boundary of the computational domain which re- 
produces most faithfully the influence of the convection zone 
on the radiative zone, while those near the lower boundary 
are selected in such a way as to impose the required dipolar 
field, but otherwise be as inconspicuous as possible. 

The set of coupled partial differential equations de- 
scribed in @ and expanded in Appendix B represents a 
12th order system, and therefore requires 12 independent 

2 although in work of Kitchatinov & Riidiger (2006), the effect of 
the flows does influence the boundary conditions applied to the 
magnetic field. 
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boundary conditions. 

Bottom boundary conditions. The bottom boundary of 
the computational domain is located at ri n , and should 
be thought of as the interface between the modelled fluid 
and an electrically and thermally conducting solid sphere 
(hereafter "the core"). Impermeability implies 



u r = at r = r- 1Xi 



(10) 



We generally impose no-slip boundary conditions at n n (ex- 
cept when specifically mentioned), so that ue = and 
u</> = fin sin 8Cli n , where f2j n is the constant angular velocity 
of the core in the rotating frame. The value of J7i n is selected 
in such a way as to guarantee that the total torque applied 
to the core be zero: 

r-rr/2 



/71 



72 



3^1 B B 

pvrfn sin 2 8 — h n n sin 8 — - — - ) sin 8d8 — 

or 4-7T 



(11) 

Assuming that the core is a homogeneous thermally 
conducting solid implies that the temperature fluctuations 
within satisfy Laplace's equation 

V 2 f in = . (12) 

The regular solution can be expanded upon standard eigen- 
functions, namely 



r n M) = ^CP n (cos0)r™ 



(13) 



where P n is the n-th order Legendre polynomial, and the set 
{a 1 ^} are integration constants. Requiring the simultaneous 
continuity of T and its derivative at r = n n while satisfying 
equation (|13|) constrains the {ajf } and yields a unique rela- 
tionship between T and dT/dr, which is used as a boundary 
condition on temperature. 

Two boundary conditions for the magnetic field are re- 
quired. They are obtained in a similar fashion to that of 
the temperature fluctuations, by assuming that the core is a 
homogeneous electrically conducting solid. In that case the 
magnetic field within satisfies 

V 2 B in = , 



which, by axial symmetry, is equivalent to 

(V 2 B in )^ = and (VxB in ) = 
Defining a flux function \ as 



Bp" = Vx 



X e 



(14) 



(15) 



(16) 



implies that both BJJ 1 and x' n satisfy the same equation, 
namely 



V 2 Si D 



— , similarly for x' n 



The regular solution for can be expanded as 
4"M) =J2bn^(Pn(oo S 8))r n , 



(17) 



(18) 



where as before the {6Jf } are integration constants. Requir- 
ing the simultaneous continuity of B$, and its derivative at 
r = rin while satisfying equation (|18[) yields a unique rela- 
tionship between and dB^/dr, which is used as one of 
the two required boundary conditions on the magnetic field. 



A point-dipole located at r = is required to main- 
tain the primordial magnetic field in this steady-state study. 
Note that this implicitly assumes that the dynamics of the 
tachocline occur on timescales much shorter than the global 
magnetic diffusion time (which is a reasonable assumption 
since the global magnetic diffusion timescale is of the order 
of the age of the Sun). To match B'p to a point dipole of 
amplitude Bo at r = r ln and 8 = (on the polar axis), we 
require that 



in D ffnSin< 
X =B — 7T 



08 



(19) 



Requiring the simultaneous continuity of B r and Be at 
r = rin provides a relationship between these two quantities, 
which is used as the second of the two required boundary 
conditions on the magnetic field. 

Top boundary conditions. The choice of boundary condi- 
tions near the top boundary should model the effects of 
the convective zone dynamics on the radiative-convective 
interface. 

Requiring the continuity of the angular velocity pro- 
file across the interface is a natural choice for the boundary 
condition on u^, although one must bear in mind that it 
is intrinsically equivalent to assuming that the interface is a 
no-slip boundary. For consistency, a similar condition should 
then also be imposed on the other component of the merid- 
ional flow parallel to the boundary, ug. Hence we select 



ue{r ou t,8) = u° e ut (8) , 

u<t>(r ou t,0) = r ou t sin 8 (f2 cz (8) - O) 



(20) 



where UQ Ut (8) could be any desired latitudinal velocity pro- 
file, and fi cz is given by equation {T}. 

Boundary conditions are also needed on the vertical ve- 
locity u r . In previous models of the radiative zone (Paper 
I; BZ06), the radiative-convective interface was assumed to 
be impermeable. However, this boundary condition is inade- 
quate as a model of a (permeable) fluid interface (see Section 
2.2|l . Instead, we consider the more general distribution: 



M r (Tout 



9) 



(21) 



where the profile selected for u° ut (8) is discussed in more 
detail in Section \4. 2. II 

The temperature boundary condition near the upper 
boundary is selected to model a very efficiently conducting 
fluid, by requiring that 



(22) 



outward of the computational domain. Matching T to T° ut 
provides a unique relationship between T and dT/dr at the 
outer boundary. 

The boundary conditions for the magnetic field are ob- 
tained by assuming that the convection zone is an infinitely 
diffusive medium (so that V 2 B° ut = outward of r out ) and 
selecting the solution which decays as r — > oo: 



B$* (r,8) = J2 C* (Pn (cos 8)) r 
n=l 

X out (r,0) = ]>>r^(Pn(cos0))r 



-(n+l) 



-(«■+!) 



(23) 
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Matching these solutions at r = r out with the solution for B 
within the domain provides relationships between and 
its derivative, and between B r and Bg respectively. 



3.3 Non-dimensional parameters, and the nature 
of boundary layers 

3.3.1 Enhanced diffusivities 

The typical solutions of the set of equations and boundary 
conditions described in Sections 13.11 and 13.21 are known to 
contain thin boundary layers and internal layers of various 
kinds, typically scaling as some (positive) power of the dif- 
fusivities. The limited achievable numerical resolution of the 
algorithm used implies that the set of equations <(9j cannot 
be solved for solar values of these quantities. To address the 
problem we multiply the diffusivities ~P(r), rj{r) and 



(24) 



by the enhancement factors /„ , f v and f K respectively. When 
these factors are all selected equal to each other the Prandtl 
and magnetic Prandlt numbers (Pr = V/k and Pr m = V/rj 
respectively) are solar everywhere. Note that we specifically 
do not view the numerical model diffusivities as "turbulent 
diffusivities" . We strive instead to derive characteristic scal- 
ings of the solutions with each of the /—factors, and under- 
stand the asymptotic behaviour of the solutions as they are 
simultaneously reduced towards unity. 

3.3.2 Non-dimensional parameters 

A better understanding of the scaling behaviour of the nu- 
merical solutions as /„ , /,, and / K are varied can be achieved 
by working with non-dimensional parameters. In what fol- 
lows, we scale all distances with the solar radius 7?©, and 
all velocities with _Rof2 eq (see Appendix C for details). The 
magnetic field strength is scaled by Bo, which is the am- 
plitude of the imposed field on the inner boundary on the 
polar axis (see equation US})- Within this framework, the 
following non-dimensional numbers naturally emerge: 



E v 



E v = 



fvV f P 



E K 



Pr 



H = 



Rq ^cq 



= UEl 







E K /« ' 



Bu = 



E V E V 
A = 

ND P 



//^H where A : 



B 2 



Bu 1 ' 



(25) 



where D p is the local density scaleheight (D p = 8.6 x 10 cm 
in the tachocline region) and iV is the buoyancy frequency 
(N = 8 X 1CT 4 s _1 in the tachocline region). All of the 
numbers denned above actually vary with position within 
the solar interior. The first three are Ekman numbers, de- 
fined as the ratios of the rotation timescale to the respective 



diffusion timescales. The solar Ekman numbers calculated 
using the microscopic diffusivities at the base of the con- 



vection zone are E 







2 x 10" 



E% = 3 x 10" 



and 



E® = 10 -9 respectively, all well below unity. The fourth is 
the Prandtl number, equal to 2 x 10 -6 near the base of the 
convection zone in the Sun. The Hartmann number H de- 
pends naturally on the magnetic field strength B selected. 
For B — 1G, H = 5 x 10" 9 at the base of the convec- 
tion zone, and H scales as the inverse of the magnetic field 
strength. The Burger number Bu is of the order of 10 3 . 



3.3.3 Expected boundary layers and desirable parameter 
range 

To run simulations in a parameter regime that is at least 
in the same asymptotic limit as the Sun, one should bear 
in mind the following constraints. To be in the non-diffusive 
regime, one must ensure that all three Ekman numbers in 
the simulations are well-below unity. The /—parameters se- 
lected for the numerical calculations should therefore remain 
below (_E ) -1 , (_E®) _1 and (_E®) _1 respectively, and also 
ensure that the hierarchy E v <C E v <C E K is respected. 

In the non-magnetic case, GB08 showed that the char- 
acteristic lengthscales of the dynamics of the radiative zone 
(for solar parameter values) can be summarised as l\ ~ Rq 

(i.e. the entire interior), I2 ~ Dp (i.e. h ~ 3 



j^Rq) and 



1/2 

h ~ E v Rq. Hence, to retain the same ordering of length- 
scales as in the Sun (I2 > h h), the factor f K / '/„ must 
not drop below 1/3, placing very strong constraints on the 
Prandtl number used in the simulations (i.e. Pr should not 
exceed ~10Pr Q ). 

The typical magnetic boundary layers are reviewed by 
Acheson & Hide (1973). In the case of a magnetic field par- 
allel to the boundary considered, the nature of the boundary 
layer depends on the respective sizes of E v and H: if E v > H 
then the boundary layer thickness is of the order of H 1 ' 2 ^©, 
while if E v < H then the boundary layer thickness scales as 
E^ 2 Rq (i.e. the boundary layer is actually of Ekman type). 
We conclude that for realistic field strengths, the bound- 
ary layers in the horizontal field case (e.g. near the equator) 
will always be of Ekman type. When the magnetic field is 
oblique, the situation changes: the Ekman regime is recov- 
ered when E v <<iC H 2 , while the boundary layer has a mag- 
netic nature when > H 2 , in which case its thickness is 
of the order of H_Rq . We therefore conclude that unless the 
magnetic field strength (in the tachocline) has amplitudes 
lower than a mG, one may actually expect a Hartmann layer 
near the radiative-convective interface. 

To guarantee that Z?„/H 2 3> 1 in the simulations we 
need f v < £ /(H ) 2 ~ 10 2 if the magnetic field is indeed 
of the order of 1G as suggested by GM98. This is clearly 
not an achievable goal, so we must inflate the magnetic field 
strength artificially to be of the order of IT near the outer 
boundary; in that case the constraint on /,, is much less 
stringent, with /,, ^ 10 10 . 

Finally, if we rely on flows downwelling from the convec- 
tion zone to confine the magnetic field, then the amplitude 
of the flows assumed at the radiative-convective interface 
must be large enough for their magnetic Reynolds number 
R m to be larger than unity. Since R m is very loosely defined 
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as 



where \u r \ is the typical amplitude of the radial flows and 
A ~ O.O3i?0 is the observed thickness of the tachocline 
(both in cgs units), then equation (|26[) implies that we 
need a priori flows with amplitude of the order of f v rj/A ~ 
2000(/„/10 10 )cm/s to have a magnetic Reynolds number of 
order one. For values of ~ 10 10 (typically achieved in the 
simulations), this naturally requires unphysically large flow 
amplitudes. The consequences of this choice are discussed in 
Sections 14.31 and [5] 

3.4 Numerical method 

3.4-1 Method and tests. 

The numerical method selected for the solution of this el- 
liptic system is based on the expansion of the governing 
equations and boundary conditions upon a truncated set 
of orthogonal functions spanning the direction (Chebi- 
shev polynomials to be precise), followed by the solution of 
the associated set of ordinary differential equations (ODEs) 
in the r-direction using a second-order Newton-Raphson- 
Kantorovich (NRK) relaxation algorithm. The numerical al- 
gorithm is globally similar to the one used by Garaud (2001) 
and in Paper I, and is summarised in Appendix C. The high 
latitudinal resolution required to capture the detailed struc- 
ture of the interior dynamics demands a high-order expan- 
sion in 8. This in turn requires a sophisticated parallelized 
version of NRK An outline of RELAX - the parallel NRK 
algorithm freely available upon request to P. Garaud - is 
given in Appendix D (see also Garaud & Garaud, 2007). 

The complexity of the full set of equations © precludes 
the derivation of analytical solutions in the general case, and 
hence direct tests of the full numerical algorithm. However, 
we have extensively and successfully tested the code through 
direct comparisons of the numerical solutions with analytical 
solutions of selected subsets of the system: Paper I tested 
the un-stratified magnetic version of the code, while GB08 
tested the stratified non-magnetic version of the code. 

3.4-2 Note on the convergence of the solutions 

In the NRK/RELAX algorithm, a guess is required for the 
solution of the ODE system, and refined by successive itera- 
tions until the desired accuracy is reached. The convergence 
of the algorithm is guaranteed and immediate for linear sys- 
tems regardless of the initial guess, but requires more accu- 
rate initial guesses for increasingly dominant nonlinearities. 
As a result, solutions are found in practise, first for high- 
values of the diffusivities (i.e. high /) and with a low num- 
ber of latitudinal modes, and then progressively followed 
into the nonlinear regime by lowering / and increasing the 
number of modes. 

Convergence of the solutions usually becomes painfully 
difficult as / is lowered below a certain threshold. The causes 
of the problem can be varied: in some cases, the nonlinear- 
ities begin to dominate the solution and bifurcations occur 
which reduce the basin of attraction of the solution followed; 
in others, the spatial structure of the solution becomes finer 



and finer, requiring more latitudinal modes and more radial 
meshpoints to be fully resolved. 

The only way to address the first problem is to follow 
the solution carefully, and reduce / only by a very small 
fraction between each run. This has the disadvantage of be- 
ing a frustratingly slow process, with little to be gained in 
practise. As a matter of fact, the likely change of stability 
of the solution tracked (suggested by the change in conver- 
gence properties) is an interesting dynamical information in 
itself (see Section f3. 4. 3|l . 

To address the second problem typically requires in- 
creasing the number of modes and meshpoints, which leads 
to a very rapid increase in computational costs. Nonethe- 
less, some acceptable trade-offs are found with typical sim- 
ulations using a few thousand meshpoints (mostly concen- 
trated in the boundary layers) and about 60-100 latitudinal 
modes. These require no more than a few hours per iteration 
on about 2 7 -2 8 processors (see Appendix D). 

3-4-3 Discussion of the selection of the algorithm 

At this point it may be worth discussing an oft-raised ques- 
tion about the selection of this specific numerical method of 
solution, as an alternative to the time-dependent ASH code 
used by BZ06. Both methods, when taken individually, are 
equally useful tools for studying the problem, each with their 
own advantages and limitations. 

For example, the algorithm of BZ06 studies the general 
properties of the radiative zone under the forcing consid- 
ered, and in addition provides information on the stability 
of the system, naturally taking into account transport by 
smaller-scale three-dimensional flows when necessary. How- 
ever, the very lengthy integration time required implies that 
only very few different parameter values can realistically be 
explored. Moreover, the time-dependent nature of the algo- 
rithm makes it more difficult to identify quasi-steady states 
(in particular when the typical timescales in the numerical 
model are far less obviously separated than in the Sun). 

By contrast, the simplified axisymmetric and steady- 
state nature of the equations solved in Paper I and here 
implies that each solution can be found in at most a day of 
real-time computations (more typically an hour). As a re- 
sult, we are able to explore a wide range of parameters and 
study in detail the behaviour of the solutions as a function 
of the imposed forcing (boundary flows, differential rotation, 
magnetic field strength, ...) and as a function of the back- 
ground model (diffusivities, thermo-dynamical profile, ...). 
This feature is invaluable when trying to understand how so- 
lutions scale with the various parameters - see GB08 for ex- 
ample. In addition, we always know that the solutions found 
are indeed the steady-state solutions of the system, and not 
transients. On the other hand, taking into account trans- 
port by three-dimensional turbulence is not possible (with- 
out artificial parameterization of its effects). In addition, we 
cannot a priori know whether the steady-states found are 
stable, or whether they may be unstable to either 2D or 3D 
perturbations. 

Interestingly, however, we can make use of the fact 
that the emergence of unstable solutions implies a reduc- 
tion of the basin of attraction of the stable ones (see Section 
13 .4. 2 p . and (tentatively) identify sudden convergence diffi- 
culties with the presence of a bifurcation. For example, the 
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simulations of BZ06 can be compared with ours for 

/„ = 2.6 x 10 s , /„ = 8 x 10 7 and f K = 8 x 10 5 , (27) 

(see Section I4.1|l . Our numerical algorithm has little diffi- 
culty converging to a solution for parameters /„, f v and 
f K respectively about 2 times larger than the ones quoted 
above. To go from there down to the values used by BZ06 on 
the other hand requires much more care, a problem which is 
very likely related to the fact that the BZ06 solution is seen 
to be mildly unstable to 2D and 3D perturbations. 



4 NUMERICAL EXPERIMENTS 

The versatility of the numerical algorithm constructed com- 
bined with its rapid real-time execution (see Appendix D) 
enables us to perform an extensive range of simulations, both 
in terms of the parameter regime studied and in terms of the 
forcing applied. In this section we use this tool in a number 
numerical experiments. 

In Section 14.11 we first prove our diagnostic for the 
failure of previous numerical attempts (Paper I, BZ06) to 
reach a dynamical equilibrium qualitatively similar to the 
one studied by GM98. We reproduce a quasi-steady state 
similar to the one achieved by BZ06, and study its scal- 
ing properties. We show that the meridional flow velocities 
induced in the tachocline, when the radiative-convective in- 
terface is modelled as an impermeable boundary, are always 
too low to confine the field. 

We contrast this result in Section 14.21 with a similar 
calculation in which meridional flows are directly pumped 
through the boundary (through some unspecified radial forc- 
ing mechanism) . We show that their amplitude remains suf- 
ficiently high in the tachocline to provide significant field 
confinement. 

We then perform a suit of simulations in Section 14.31 
in the "desirable" parameter regime specifically selected in 
Section f3. 3. 31 and discuss the numerically derived solutions 
both qualitatively and quantitatively. We show that the sys- 
tem can, in some situations and for low-enough diffusivities, 
achieve a balance qualitatively similar to that of the GM98 
model, but that the solutions also reveal a number of pre- 
viously unaccounted-for dynamical subtelties (see Sections 
IP | KE\ and [K3 > . 



4.1 Failure of "impermeable" simulations 

4- 1.1 The calculation 

In this first numerical experiment, we compare our quasi- 
steady solutions with those integrated by BZ06. Note that 
since the overall interior magnetic field strength decays 
slightly in their time-dependent calculations, while our 
steady-state solutions fix the value of the field amplitude 
near the inner boundary, the two simulations can only be 
compared qualitatively at best. Nonetheless, the compari- 
son is meaningful since the system is thought to relax to a 
quasi-steady equilibrium on a timescale comparable with the 
meridional circulation timescale, which is argued (by BZ06) 
to be much shorter than the magnetic diffusion timescale. 




Figure 2. Quasi-steady state results for a numerical simulation 
similar to that of Brun & Zahn (see Section [4.11 . In each quadrant, 
the dotted circle marks the position of the base of the convection 
zone at r = r cz , the outer solid circle marks the solar surface 
and the inner solid circle is at r = r ln . Top left: Streamlines of 
the meridional flow. Clockwise flows are shown as solid lines, and 
anticlockwise flows are shown as dotted lines. Top right: Toroidal 
field (in Gauss). Bottom left: Magnetic field lines. Bottom right: 
Angular velocity profile. Note how the magnetic field lines connect 
to the bottom of the convection zone and establish a differentially 
rotating Ferraro state throughout the domain. 

To achieve the same parameter regime as the one se- 
lected by BZ06 in terms of the diffusivities, we set 

/„ = 2.6 x 10 s / , /„ = 8 X 10 7 / and / K = 8 x 10 5 / (28) 

and progressively reduce / from about 10 5 to 1. When / = 1, 
our diffusivities are the same in the tachocline region as 
their^EI- Naturally, we select u° ut = to mimic an imper- 
meable boundary as in BZ06. We also modify our bottom 
boundary conditions to be "stress-free" , again to be consis- 
tent with their simulations. Finally, the radial component of 
the magnetic field on the inner boundary is set to be same 
as the initial condition selected by BZ06 (i.e. B r — 500G at 
the poles at r = n n ). 

4-1-2 Qualitative comparison of the results 

The results for / = 1 (see equation (|28jl ) are shown in Figure 
[2] The overall structure of the streamlines and of the mag- 
netic fields lines, as well as the amplitude of the toroidal 

3 Also note that the BZ06 simulations actually have constant dif- 
fusivities throughout the radiative interior, while ours vary with 
radius as in the Sun. Since we only seek a qualitative comparison 
of the two simulations, this difference in the diffusivities near the 
inner boundary is irrelevant. 
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As a result, one can straightforwardly see that the mag- 
netic Reynolds number does not increase with decreasing 
E v when the respective diffusivity ratios are held constant: 
with R m defined as in equation (|26|l , if \u r \ is proportional 
to V f K fv then the magnetic Reynolds number is in fact pro- 
portional to \/f v f K /fr} and remains roughly constant as / 
decreases. The effective magnetic Reynolds number of the 
BZ06 simulations can be estimated from the flow amplitudes 
calculated in Figure0to be R m ~ 5 x 10 -3 - far too low to 
influence the magnetic field. 

Using the true solar diffusivities would not change the 
situation since in that case it can be shown that the expected 
magnetic Reynolds number is of the order of R m ~ 4x 10~ 2 . 
We have therefore confirmed our original diagnostic, namely 
that the reason for the failure of numerical simulations to 
find confined-field solutions is not related to the fact that the 
diffusivities in the simulations are too high (i.e. selecting so- 
lar diffusivities would yield a similar result), but rather to a 
fundamental mis-representation of the radiative-convective 
interface. 



Figure 3. Absolute value of the radial component of the merid- 
ional flow at r = O.68i?0 as a function of latitude, when / = 1 
(plus symbols), / = 10 (diamonds), and / = 100 (triangles), for 
/ defined in equation {28}. The radial position was selected to be 
safely below the Ekman layer for all three calculations. The solid 
part of the curve denotes u r > and the dotted part of the curve 
denotes u r < 0. Note how \u r \ roughly scales with /. 

field and of the differential rotation compare qualitatively 
well with the simulations of BZ06 (see the t = 4.7Gyr solu- 
tion in Figures 4 and 5 of BZ06) . Small discrepancies in the 
magnetic field geometry can be attributed to the fact that 
their magnetic field is diffusing out of their inner core, while 
it is maintained in our simulation by a permanent dipole. 

This steady-state calculation therefore confirms the 
main conclusion of the study carried out by BZ06, namely 
that in the case where the radiative-convective interface is 
modelled as an impermeable boundary, the only possible 
quasi-steady solution has an unconfined field structure lead- 
ing to a differentially rotating Ferraro state. 

4-1.3 The nature and amplitude of the flows 

We now justify our diagnostic of the failure of "impermeable 
simulations" to find a confined-field solution. The induced 
meridional flows in this set of simulations are generated pri- 
marily by Ekman pumping on the boundary. Indeed, the 
magnetic field strength selected by BZ06 implies a field of 
the order of a few tens of Gauss at the most near the outer 
boundary, and therefore according to the analysis of Section 
13. 3. 31 a boundary layer in the Ekman regime rather than in 
the Hartmann regime. In that case, the flow amplitudes are 
in fact well-estimated by the (non-magnetic) work of GB08 
and scale as _E„/\/Pr. This is verified in Figure [3] which 
shows the radial velocity profile of the flows as a function 
of latitude, well below the Ekman layer (at r — O.68-R0), in 
three different simulations: for / = 1, / = 10 and / = 100. 
It is clear from the figure that the meridional flow velocity 
decreases linearly with /, and can in fact be shown to scale 
as predicted above. 



4.2 The possibility of magnetic field confinement 
by the t her mo- viscous mode 

Following the original idea proposed by GM98, we now ex- 
plore the possibility of confining the interior field through 
flows which are directly downwelling from the convection 
zone. As shown by GB08, flows which are pumped into and 
out of the radiative zone (by stresses in the convection zone 
for example) can retain significant amplitudes throughout 
much of the tachocline (and below), and could therefore re- 
sult in a magnetic Reynolds number larger than unity for low 
enough diffusivities. We illustrate this idea with a qualitative 
example, which, for ease of comparison with the results of 
the previous section, is selected to have the same diffusivity 
ratios and magnetic field strength as those of BZ06. 



4-2.1 Structure of the imposed meridional flow 

We now select a radial velocity profile u° ut (0) which satisfies 
the following properties: (i) u° ut (0) = (the radial velocity 
is zero on the polar axis), and the flow is symmetric with 
respect to the equator (ii) the total mass flux through the 
boundary (at r — r cz ) is zero: 

«r/2 

/ u° ut {8)sm9d8 = , (29) 
Jo 

and (iii) the total angular momentum flux carried by the 
meridional circulation through the boundary is zero: 

rn/2 

/ r 2 sin 2 8D. cz (8)u° ut (9)d8 = . (30) 
Jo 

The simplest solution which satisfies all three constraints is 
a 6-th order polynomial in /1 = cos(#): 

uTM = U (l + 53 a nf i 2n ^j (31) 
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Figure 4. Example of imposed radial flow u° ut as a function of 
latitude, for the values of a 2 and 04 described in equation ((TJ, 
satisfying the conditions listed in the main text. The velocity of 
the flow (Uo) here is arbitrarily selected to be —1 cm/s at the 
equator. The width of the upwelling region, as can be seen in the 
figure, spans about 30° in latitude and is centred on 30° . 



where 



cti = 



a.2 



'<:■ = 



3(715 + 169a 2 + UI04) 
143 + 13a 2 + 25a 4 
5(1001 + 299a 2 + 207a 4 ) 
143 + 13a 2 + 25a 4 
91(33 + lla 2 + 7o 4 ) 
143 + 13a 2 + 25a 4 



(32) 



and Uo is the radial velocity at the equator. 

An example of this profile can be seen in Figure |4] for 
the values of a 2 and 0,4 given in ^ the upwelling region is 
centered at 30° latitude, and ranges from 15° to 45°. Note 
that the selection of the "simplest" profile is fairly arbitrary 
- had we chosen a different set of orthogonal functions, we 
would have obtained a slightly different flow profile. One 
may also wish instead to select a profile where the width 
and position of the upwelling region is prescribed, or in- 
stead of selecting the polar velocity to be zero, to select 
the equatorial velocity to be zero, etc... The only necessary 
constraints, however, are zero mass and angular momentum 
fluxes. These conditions are often ignored, but are crucial to 
obtaining meaningful results; failure to satisfy them yields 
unphysical angular velocity profiles (see Garaud, 2007). 



4-2.2 Qualitative structure of the solutions 

In the simulation presented here, the diffusivities are se- 
lected as in equation (|28p with / = 1. A radial meridional 
flow profile is imposed at the outer boundary (see equation 
(|31[) ). with an amplitude Uo which is progressively increased 
from to about -5.5 cm/s. The resulting numerical solution 
is presented in Figure [5] It clearly shows that the field lines 
are strongly distorted by the meridional flows, and seem to 
be confined both near the equator, and more importantly 
also near the poles. In fact, our results appear to be at least 
qualitatively similar to the analytical solutions of Wood & 
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Figure 5. Same as [2] but with an imposed meridional flow at 
r = r out (see Figure g}. 



Mclntyre (2007), which focussed on the polar regions only. 
The angular velocity profile is very different from the one 
observed in the simulations shown in Figure [2] and reveals 
uniform rotation from the equator up to about 60° in lat- 
itude. The rotation rate of the inner core is found to be 
0.86n oq . 

We have therefore shown that, at least on a qualitative 
basis, magnetic field confinement is possible provided a ra- 
dial flow is imposed near the radiative-convective interface. 



4.3 Exploration of parameter space 

Having proved on a qualitative basis that the GM98 model 
may indeed lead to global field confinement, we now turn 
to a more quantitative analysis of the resulting tachocline 
dynamics. For this purpose, we return to using the boundary 
conditions discussed in Section \3. 2 1 

In what will be referred to from here on as the "fiducial 
model", the enhancements factors /„, f n and f K are now 
selected as follows to be in the "desirable" parameter range 
discussed in Section [3.3.31 



f ~ f 

f»= 10 ' = /./*= 100 



(33) 



and the factor / is progressively reduced from 10 down to 
10 10 typically (8 x 10 9 for the lowest case). For / = 10 10 , 
and in the tachocline region, £„~2x 10 -6 , £,~3x 10~ 4 
and Ek ~ 10 _1 . In that case, and for all simulations, 
Pr = 10Pr s (guaranteeing that the thermo- viscous mode 
lengthscale is indeed of the order of the entire radiative 
zone), and Pr m = OTPr^. Note that the simulations of BZ06 
use a Prandtl number of the order of 300Pr®, which does 
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Figure 6. Evolution of the poloidal field topology and of the an- 
gular velocity profile in the fiducial model for, from top to bottom, 
/ = 10 13 , 10 12 , 10 11 and 10 10 . Note that regions with Q < 0.6Q eq 
are drawn in black. 



not satisfy our criterion for the hierarchy of the expected 
characteristic lengthscales. 

The magnetic field strength on the inner boundary at 
the pole is selected to be Bo = 7T, so that the field strength 
in the tachocline would be of the order of IT in the absence 
of meridional flows. The radial flow velocity at the outer 
boundary is the one discussed in Section 14.2.11 with Uo = 
— 550cm/s, and the latitudinal flow velocity is everywhere 
zero. 



Streamlines T (K) 




Figure 7. Evolution of the topology of the streamlines and of the 
temperature perturbations in the fiducial model for, from top to 
bottom, / = 10 13 , 10 12 , 10 11 and 10 10 . Note that the streamlines 
shown are representative streamlines, the contours selected are 
different in each plot. Also note that the temperature scale is 
different in each plot. 



4-3.1 From the diffusive regime toward the asymptotic 
regime 

The qualitative evolution of the poloidal field topology and 
of the angular velocity profile as / is reduced can be seen 
in Figure [6] while the corresponding figure for the stream- 
lines and the temperature profile is Figure [7] In all simu- 
lations, the viscous Ekman number is much smaller than 
one. The magnetic field begins to affect the angular ve- 
locity profile as / decreases from 10 13 to 10 12 , notably in 
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Figure 8. Radial flow velocities at r = O.68B0 as a function of 
latitude in the fiducial model, for three values of /: / = 10 10 (plus 
symbol), / = 3 X 10 10 (diamonds) and / = 10 11 (triangles). As 
in Figure [3] the solid lines denote upward flow (u r > 0) and the 
dotted lines denote downward flow (u r < 0). The figure therefore 
shows a downwelling in mid-latitude, with upwelling at the poles 
and the equator. Note how the amplitude of the radial velocity 
appears to increase with decreasing /. 

the deeper interioiQ. This transition occurs when the mag- 
netic Ekman number drops below one. The structure of the 
temperature perturbation profile changes as / is decreased 
from 10 12 down to 10 11 , which is attributed this time to the 
"thermal" Ekman number approaching one. Finally, there 
is a clear transition between the "unconfined field" config- 
uration for / = 10 11 and the "confined field" configuration 
for / = 10 10 , which corresponds to the magnetic Reynolds 
number increasing to values above unity. 

Contrary to the impermeable-boundary case studied by 
BZ06, the radial flow velocities do not decrease with the 
diffusivities, but are actually seen to increase as / is re- 
duced. This is illustrated in Figure [8] which shows the radial 
flow velocity near the base of the presumed tachocline (at 
r = O.68-R0, as in Figure [3]), as a function of latitude, for 
three values of /. Provided E K < 1, it appears that u r oc 1// 
although this scaling cannot, for numerical reasons, be ex- 
plored over many orders of magnitude. As a result, the non- 
linear interaction of the flow and the field rapidly becomes 
stronger and the field appears to be more and more confined, 
as can be seen in the last panel of Figure [U 

4-3.2 The "lowest diffusivities" simulation 

We now decrease / as far as possible, until convergence be- 
comes too difficult. The lowest value achieved in the fiducial 
model is / = 8 x 10 9 and the results are presented in Fig- 
ure [5] The number of radial meshpoints used is 3000, and 

4 Recall that the magnetic diffusivity is smaller in the deeper 
interior, see Appendix A. 



Dynamics of the solar tachocline 13 



Streamlines Temperature pert. (K) 





20 40 60 80 80 60 40 20 



Figure 9. Results of the simulations of the fiducial model for 
/ = 8 X 10 9 . For each quantity plotted (from top left to bot- 
tom right, respectively, representative streamlines, temperature 
perturbations, field lines and rotation rate), we show both a full 
quadrant as well as a zoomed-in strip of outer boundary layer (the 
region between r = O.65i?0 and r = r ou t). As usual, the solid 
streamlines denote clock- wise flows, and dotted streamlines de- 
note counter-clockwise flows. Note that regions with Q < 0.6f2 C q 
are drawn in black. 



the number of latitudinal modes is 80 (note that since only 
even modes are selected to guarantee equatorial symmetry, 
solving for 80 modes means that the highest Fourier mode 
considered actually is cos(158#), see Appendix C). 

The top- right quadrant of Figure[9]shows representative 
streamlines. The Ekman layer near the outer boundary, of 
width ~ O.OOli?0, is just visible in the top-right strip. The 
Ekman layer flow does indeed downwell at the poles and 
equator, and upwells in mid-latitudes. 

More clearly visible is the "thermo- viscous" mode below 
the Ekman layer, which (rather unexpectedly, to be honest) 
has the opposite structure: downwelling in mid-latitudes and 
upwelling near the poles and equator. This mode appears to 
be confined roughly within r £ [0.67, 0.7]Rq, an effect which 
can only be attributed to the Lorentz forces arising from the 
magnetic field underneath. Indeed, the solutions of GB08 in 
the absence of magnetic fields in the same parameter regime 
do not reveal any structure on this typical lengthscale. The 
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associated effect of the flow in confining the magnetic field 
is also obvious in the bottom-left quadrant: the field lines 
are most strongly distorted away from the strictly dipolar 
structure in the same region (r £ [0.67, 0.7]Rq). Field con- 
finement is discussed in more detail in Section \4. 41 

The solution for the angular velocity profile is again 
more complex than expected: the combination of the two 
modes of propagation of the meridional flows creates radial 
structures both on the Ekman scale and on a larger scale 
(r ~ O.IRq), as well as latitudinal structures which are 
clearly seen to be correlated with the magnetic field lines. 
However, contrary to the solution shown in Figure[5] there is 
no clear evidence for a large-scale latitudinal or radial gradi- 
ent in Q below about O.6i?0, including in the polar regions. 
This is the first set of self-consistent simulations to reveal 
this kind of solution. Since the angular velocity profile in- 
ferred from helioseismic inversions can be thought of as a 
weighted spatial average of the true angular velocity pro- 
file in the Sun, and since the spatial extent of the averaging 
kernels is much larger than the features seen in the simula- 
tions, the angular velocity profile shown in Figure [5] would 
be observationally interpreted to be uniform below O.6i?0, 
whereas the one shown in Figure [2] would not. The actual 
value of the rotation rate calculated in the simulations is 
discussed in Section [45] 

A closer look at the angular velocity profile in the re- 
gion r € [0.65, O.7]i?0 (see the bottom-right strip) reveals a 
change in the dominant angular momentum transport pro- 
cesses between the uppermost layers (r > 0.67-Rq) and the 
lower layers (r < 0.67Rq). When r > 0.67-R© the angular 
velocity contours are more-or-less aligned with the stream- 
lines, whereas for r < 0.67i?o the angular velocity contours 
are aligned with the magnetic field lines. This is at least 
qualitatively consistent with the idea proposed by GM98 of 
a multi-layered tachocline, in which the dynamics of the up- 
permost layer are controlled by the tachocline flows, while 
the dynamics of the bulk of the radiative zone are controlled 
by the primordial field. The interface between the two re- 
gions, which could be identified with the "magnetic diffu- 
sion layer" of the GM98 model, clearly has quite a complex 
topology. 

One can also observe a number of very strong localised 
"jets" of more rapidly or more slowly rotating fluid. The jets 
appear to be stronger for lower diffusivities (one can easily 
compare the results of Figure [9] with the last panel of Figure 
O. In fact, we tentatively attribute the numerical scheme's 
failure to find solution for / < 8 x 10 9 to an intrinsic instabil- 
ity of the jets. But whether such strong features would exist 
in the Sun is debatable: both the imposed flow velocities and 
the magnetic field strength have been artificially increased 
by many orders of magnitude to ensure a magnetic Reynolds 
number greater than one in the simulations, and a Hartmann 
number closer to that of the Sun (see Section T3. 3. 3 [) . As a re- 
sult, the amplitude of these jets is also artificially increased 
in the simulations compared with what one may expect in 
the Sun. Unfortunately, given the nonlinear nature of the 
problem, it is difficult to predict what the actual strength of 
these jets would be should diffusivities, magnetic field and 
imposed flows velocities be simultaneously reduced to the 
expected solar values. 



4.4 Magnetic field confinement 

To study the nature of field confinement more quantitatively, 
we consider the ratio £ = |-B r /_B r:d ipoiar| where -B r ,di P oiar is 
the hypothetical radial component of the solution for the 
magnetic field in the purely diffusive case: 

B r , dipolar (r, 9) = Bo cos 9 (j-^ . (34) 

Typically, we expect that £ < 1 when magnetic field lines are 
either completely expelled from a region, or bent to the hor- 
izontal. On the other hand, we expect that £ > f in regions 
where the magnetic field lines are pushed together by con- 
verging (latitudinal) flows. It is therefore a good diagnostic 
of the processes we are interested in. 

The quantity £ is shown in Figure [TDJfor various simula- 
tions. The top two quadrants are derived from the numerical 
solutions for / = I0 11 and / = 8 X I0 9 respectively, and re- 
veal a strong change in the behaviour of the solutions when 
R m < 1 and R m > I. 

In the case where / = 10 11 , we observe a good match 
between the spatial variation of £ and the imposed forc- 
ing: £ < 1 in downwelling regions and £ > 1 in the up- 
welling regions (at mid- latitudes). However, the amplitude 
of £ remains close to one, revealing only a very weak ef- 
fect of the flow on the field. This is not entirely surprising 
since the magnetic Reynolds number for / = 10 11 can be 
deduced from the flow amplitude shown in Figure [S] to be 
Rm ~ 0.03 at r = 0.68-Rq. In fact, in this very weakly non- 
linear regime there is some evidence for a correlation (see 
Figure [TOl bottom left quadrant) between 1 — £ and /: we 
find that 10 10 /" 1 . 

However, Figure \W\ also reveals that matters become 
more complex as R m increases above unity. For / = 8 x 10 9 
the magnetic Reynolds number is now of the order of R m ~ 
1.5, and the connection between the field and the meridional 
flow is correspondingly fully nonlinear. The correlation be- 
tween 1 — £ and / breaks down, as seen in the bottom right 
quadrant. This incidentally proves that the kind of bound- 
ary conditions advocated by Kitchatinov & Rudiger (2006) 
to mimic "field confinement by meridional flows" cannot ac- 
curately describe the nonlinear dynamics of the system. We 
also observe that the regions of enhancement and reduction 
of the radial field are no longer perfectly correlated with 
the imposed upwelling/downwelling regions. In fact, some 
enhancement in the radial field strength can be seen near 
the equator as well as near the polar regions. This can be 
attributed to the convergence points of the meridional flows 
below the Ekman layer, which, as seen in Figure [§] have the 
opposite sign as that of the imposed flow. 

All in all, although there is definite evidence for mag- 
netic field confinement along the lines of the model proposed 
by GM98, it is probably fair to say that the system be- 
haves in a far more complex way than anticipated, and that 
the simulations in the fiducial model are just beginning to 
scratch the surface of the "interesting" regime. 

4.5 Interior angular velocity 

As in Paper I, we consider the angular velocity of the inner 
core ilin = ^ + fiin (which is an eigenvalue of the calculation 
performed) as a diagnostic of the most important angular 
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Figure 10. Measure of magnetic confinement, using the quantity 
£ as defined in the main text, Section 14.41 Top left: Variation of 
£ with radius and latitude, in a fiducial model for / = 10 11 . Top 
right: Same as on the left but for / = 8 X 10 9 . Note the change in 
the scale of the perturbations between the two plots, and also note 
that here the contours are logarithmically spaced. Bottom left: 
|1— £| evaluated at the outer boundary, at three different latitudes 
respectively in the downwelling regions (10°, plus symbols, and 
60°, diamonds) and in the middle of the upwelling region (30°, 
triangles). This log-log plot emphasises the power-law relationship 
between 1 — £ and / for R m < 1, see main text. Bottom right: 
Similar to the left-side plot but showing £ with a linear scale, 
emphasising the breakdown of this relationship when R m > 1. 



momentum transport processes, and of the "accuracy" of 
the model in reproducing the observations. This quantity is 
shown in Figure [TT] as a function of / for the fiducial model, 
and for simulations with faster and slower imposed flows. 

In very diffusive simulations, the predicted angular ve- 
locity f2j n ~ 0.957f2 cq is consistent with purely viscous angu- 
lar momentum transport throughout the entire interior (see 
Gough, 1985). As / is progressively reduced, the system un- 
dergoes a rather impressive bifurcation at / ~ 2.5 x 10 12 , 
which corresponds to the change of nature of the bound- 
ary layer near the inner core from an Ekman layer to an 
Ekman-Hartmann layer. When this happens, a new flow cell 
of opposite vorticity appears right against the inner core, so 
that angular momentum transport changes direction. This 
leads to a rather dramatic jump in the inner core rotation 
rate, which is unlikely to be of any relevance to the Sun, but 
deserved an explanation. Note that the critical value of / for 
which the bifurcation occurs is proportional to the imposed 
magnetic field strength Bq. 

In the case of the fiducial model, the angular velocity 
of the core then remains roughly constant (£2i n ~ (0.875 ± 
0.005)f2 cq ) as / is further reduced by two orders of magni- 
tude. This range corresponds to the regime where R m < 1, 
where the magnetic field lines essentially remain dipolar. In 
fact, the same inner core velocity is also found for lower 
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Figure 11. Angular velocity of the inner core as a function of /, 
in the fiducial model (star symbols), and in a model where the am- 
plitude of the imposed radial flow is ten times larger (triangular 
symbols) and ten times smaller (diamond symbols) respectively. 



imposed flow velocities (see Figure fTTj) . as well as for higher 
and lower field strength (always with R m < 1). This suggests 
that the value fii n ~ 0.875fi cq is a rather universal charac- 
teristic of angular momentum transport by dipolar fields in 
the Hartmann regime, in this particular geometry, and sub- 
ject to this particular imposed angular velocity profile at the 
outer boundary. 

However, it is clear from the analysis of Section 14.41 
that the simulations in the fiducial model have not quite 
yet reached the asymptotic regime even for the lowest val- 
ues of / achieved. As a matter of fact, in the last few points 
of the fiducial model curve (in the region where R m begins to 
exceed one), Q ln is seen to change more rapidly, decreasing 
slightly below 0.87 f2 cq . 

The role of the meridional flows as angular momentum 
transporters can best be seen in a simulation with much 
higher imposed flow velocities (see the curve with the trian- 
gular symbols, for Uo = — 5500cm/s). In these simulations, 
the transition to R m > 1 occurs for values of / typically 
10 times larger than in the fiducial model. In that case, we 
observe that the predicted core velocity is notably different 
from the other two cases, and changes rather rapidly with 
/. It is therefore more than likely that the predicted angular 
velocity of the core in the fiducial model would also begin to 
deviate more noticeably away from Q ln — 0.875-Rq were we 
able to continue decreasing /. Again, we are only beginning 
to approach the asymptotic regime so that the simulated 
angular velocity of the core cannot yet and should not be 
compared directly with that of the Sun. 
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5 DISCUSSION AND PROSPECTS. 
5.1 Summary of the results 

Despite the difficulties encountered when attempting to find 
numerical solutions for asymptotically low values of the dif- 
fusivities, the nonlinear dynamics which emerge from our 
simulations are closer to what one may expect from the 
GM98 model than any other simulation performed to date. 

More precisely we do observe (see Section 14. 3p the 
quenching of the large-scale differential rotation by the pri- 
mordial magnetic field, even in the polar regions. We also 
observe the partial confinement of the field by the merid- 
ional flows to the radiative interior (see Section |4.4[) . with 
the reduction of the radial field strength (in some regions) 
by more than 70%. Finally, we observe the concurrent con- 
finement of the meridional flows to the upper layers of the 
radiative zone (r > 0.67-Rq) by the magnetic field. We have 
not yet been able to reduce the diflusivities down sufficiently 
far to observe a truly segregated structure where the bulk 
of the tachocline flows are completely magnetic-free. How- 
ever, there is evidence in the observed rotation profile (see 
Figure [9} for a transition between regions where angular- 
momentum transport is dominated by the meridional flows, 
and regions where it is dominated by the magnetic field. 
The "magnetic diffusion layer" studied by GM98, which is 
thought to control this transition, appears to have a rather 
complex geometrical structure which prevents a more de- 
tailed study of the GM98 scalings. 

The calculated value of the interior angular velocity Qi n 
(see Section I4.5p does not match the observed value. How- 
ever, since even the lowest-diflfusivity simulations presented 
here are only just beginning to enter the asymptotic pa- 
rameter regime described in Section 13.3.31 we do not view 
the poor match with the observations as an intrinsic prob- 
lem with the GM98 model (yet) but rather as evidence that 
more work should be done to decrease the diflusivities even 
further - a challenging task. 



5.2 The stability of the solutions? 

As described in Section [473] the convergence of the solutions 
becomes rather difficult when the values of the diflusivi- 
ties are decreased below a certain threshold (in the fiducial 
model for / < 8 x 10 9 ). The typical reasons for this con- 
vergence failure were listed and discussed in Section 13.4.21 
intrinsic linear/nonlinear instabilities or insufficient latitudi- 
nal resolution. In this particular case, the various scenarios 
are equally plausible, and difficult to disentangle without a 
supporting fully nonlinear 3D calculation. The field configu- 
ration near the inner core could be subject to instabilities (as 
seen in the simulations of BZ06). The Ekman(-Hartmann) 
layer near the outer boundary could also be subject to insta- 
bilities, or could be developing latitudinal structures which 
are too fine to resolve (see Figure [12}. We have in fact ten- 
tatively identified the emergence of strong jets in the angu- 
lar velocity profile (see Figure [9] for example) as the most 
likely source of instabilities (and convergence failure), but 
the reader should be cautionned that this statement is only 
speculative. 

As mentioned in Section 14.31 whether the equilibrium 
governed by the GM98 model is stable or unstable cannot 
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Figure 12. Angular velocity profile on and just below the 
radiative— convective interface (a) at 0.7-Rq (b) at 0.6999-Rq, 
(c) at 0.6995-Rq and (d) at 0.699-Rq in the fiducial model for 
/ = 8 X 10 9 . Note how the angular velocity profile in (and there- 
fore also below) the Ekman layer is very different from the one 
imposed by the convection zone. 

directly be inferred from the stability of the numerical so- 
lutions (since they operate in a different parameter regime) 
but it is clearly a fundamental question. Should our diagnos- 
tic concerning the stability of the jets prove to be correct, 
then a possible path for future investigation (using this re- 
laxation method) may be to continue searching for solutions, 
starting where we left off, and slowly lowering simultane- 
ously the amplitude of the imposed flow (and field) with 
the diflusivities to limit the strength of the jets generated. 
This could be seen as blindly navigating the stable regions of 
parameter space while avoiding instability reefs. But what 
could ideally be derived from such an exercise, should we 
be able to acquire sufficient data, are scaling relationships 
between the jet strength, the diflusivities and the imposed 
flow strength which may then be used to infer tentative in- 
formation on the stability of the GM98 solution itself. 

5.3 The role of an Ekman layer? 

Our simulations have also revealed a fundamental difference 
with the original GM98 model: the presence and role of an 
Ekman mod^B- While it is clear from Figure [9] that the Ek- 
man flows themselves play no role in confining the field (con- 
trary to the claims made by Kitchatinov & Riidiger 2006) 
our simulations reveal that the role of the Ekman layer is 
still far from trivial. 

Indeed, the differential rotation profile imposed at the 
radiative-convective interface is quite different from the dif- 
ferential rotation profile transmitted by the Ekman layer. 

5 Since the bulk of the GM98 tachocline is presumed to be 
magnetic-free, one should indeed consider an Ekman mode rather 
than an Ekman-Hartmann mode 
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This is illustrated in Figure 1121 which shows the evolution 
of the angular velocity profile with depth across the Ekman 
layer. Thus while the GM98 model correctly describes the 
non- viscous tachocline dynamics below the Ekman layer, the 
Ekman mode could in fact influence the system by modify- 
ing the angular velocity profile "seen" by the bulk of the 
tachoclinqj. 

The extent to which this effect influences the tachocline 
is unclear. Firstly, other processes which also transport 
angular momentum (convective overshoot, gravity waves, 
small-scale and large-scale magnetic stresses associated with 
the dynamo field) are known to take place in the close 
vicinity of the radiative-convective interface. These pro- 
cesses were not included in GM98's analysis and cannot be 
modelled with the present numerical algorithm, but could 
equally affect the tachocline dynamics. Secondly, Ekman lay- 
ers (laminar or turbulent) must by definition be present in 
any rotating system which exhibits very rapid changes in 
the imposed stresses. However, whether the Ekman mode 
actually plays an important role in the tachocline dynamics 
depends on the relative importance of the bulk thermal-wind 
stresses and the combination of all the other rapidly vary- 
ing turbulent stresses (see Figure 1 13 p . Indeed, if the base of 
the convection zone is already essentially in thermal-wind 
balance (see Miesch, 2005 for instance), then the amplitude 
of the Ekman mode and its effect on the angular velocity 
profile will be small (GB08). On the other hand, if this is 
not the case then significant Ekman flows are expected (as 
in the simulations shown in the present work for example), 
with the aforementionned consequences. 
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Figure 13. A pictorial representation of the expected tachocline 
flows in two extreme cases. In both pictures, flows are downwelling 
from the convection zone into the radiative zone. As in the GM08 
model, the bulk of the tachocline is ventilated by the thermo- 
viscous mode, which interacts with the magnetic field in a thin 
advection-diffusion layer. On the left, and as in the simulations 
presented in this paper, the convection zone is not dominated 
by thermal-wind balance, and a significant portion of the flows 
downwelling from the convection zone rapidly return within a thin 
Ekman layer. The angular velocity profile seen by the tachocline 
differs from the one observed in the convection zone. On the right, 
a hypothetical situation where the convection zone is essentially in 
thermal-wind balance. In that case, the Ekman mode is negligible 
and the GM98 model directly applies. 



5.4 Prospects 

The present work has emphasized the necessity of flows 
downwelling from the solar convection zone as a means to 
confine the internal primordial field and guarantee the uni- 
form rotation of the radiative interior as originally proposed 
by GM98. The nature of the thermal and dynamical balance 
governing the convection zone itself therefore also controls 
the dynamics of the tachocline through the spatial variation 
and amplitude of the flows crossing the radiative-convective 
interface. Meanwhile, steady progress in modelling the con- 
vection zone has emphasized its dependence on the thermal 
stratification of the tachocline (Rempel, 2005; Miesch, Brun 
& Toomre, 2006). One can only conclude that the dynamics 
of the convective and radiative regions are intrinsically and 
fundamentally coupled through this internal layer called the 
tachocline, and that the only sensible way forward from this 
point on is the construction of whole-Sun models including 
both regions. 
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APPENDIX A: SELECTION OF THE 
BACKGROUND STATE 

The spherically symmetric background solar model selected 
for this paper is Model S of Christensen-Dalsgaard et al. 
(1991). Model S provides calculated solar data for all of the 
relevant thermodynamical and compositional quantities on 
a fixed mesh. In particular, we use the provided fields T, p, 
~p~, g, N (where N is the buoyancy frequency), c p and kr 
(where c p is the specific heat at constant pressure, and Kb. 
is the Rosseland mean opacity). In order to use this data 
for our purpose, we need to interpolate it upon our own 
selected numerical mesh, using a standard rational inter- 
polation routine (cf. Numerical Recipes). The task is del- 
icate since any "roughness" in the interpolated functions 
results in failure of convergence of the numerical algorithm, 
and the two numerical meshes are intrinsically non-uniform: 
Model S has meshpoints strongly concentrated near the so- 
lar surface, near r = and at the base of the convection 
zone r cz = 0.713-R©, while our numerical mesh has mesh- 
points strongly concentrated near the domain boundaries 
and Tout- A new interpolating routine was created which 
automatically selects a certain number of points from the 
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Figure Al. Background pressure, density and temperature (top 
left corner), buoyancy frequency squared (top right corner), mi- 
croscopic diffusivitics (bottom left corner) and ratio of diffusivities 
(bottom right corner) from Model S of Christensen-Dalsgaard et 
al. (1991). 



original mesh, performs the rational interpolation upon the 
desired mesh, and checks for the smoothness of the inter- 
polated function and of its derivative. Should the result be 
inadequate, the procedure is repeated with a different set of 
target points. Selected background quantities are shown in 
Figure lATI 

The quantities V, rj, and k are not provided by Model S, 
and must instead be calculated using the formulae derived 
by Gough (2007): 



v(r) = v c . 



n(r) 
k{r) 



+ 0.9 



-*■ cz / \ pcz 



T\ 5/ V -p 



KR 

KR,cz 

In A 



Tf \ - 3 / 2 



T 

pc p K(r) , 



Pcz J \ In A c 
In A 



In A c 



(Al) 



where the Coulomb logarithm In A and the thermal diffusiv- 
ity 7t are given by 

3/2 

A(r) = Ac 



p N ^ 



ft(r) 



Try 



KR,ci 



(A2) 



The following quantities are those at the base of the con- 
vection zone: T cz = 2.3 x 10 6 K, p cz = 0.21g/cm 3 , kr iCZ = 
19cm 2 /g, lnA cz = 2.5, v cz = 27cm 2 /s, r7 cz = 410cm 2 /s and 
Kcz = 1.4 x 10 7 cm 2 /s. 



APPENDIX B: MODEL EQUATIONS IN 
SPHERICAL COORDINATES 

The model equations described in the system (O are now 
expanded in a spherical coordinate system (r,6,<f)), with 
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u = (u r ,u g ,U4,), B = (_B r ,B 9 ,_B ) and j = (j T ,3e,3$) 
The r— component of the momentum equation: 



-psin6M2fT + ^ 



dp _d$ . „ . n 

■ a I U <P = ~n P~, 1" 3» B 4> - 3<t> B 

r sin o I or dr 



+ /„(v-n) r . 

The 8— component of the momentum equation: 
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The </;— component of the momentum equation: 



p 20 + 



(cos dug + sin 8u r 



= j r Bg - jgB r 

+ /„(V-II)* . (B3) 



where the role of f v was discussed in section 13.3.11 The di- 
vergence of the viscous stress tensor, for a stratified fluid, 
can be derived from Batchelor (1994 edition, pp. 147 and 
601). 

The mass conservation equation: 

±^(r 2 pu r ) + ^—-?-( a m8u e )=0. (B4) 
r 2 or v ' r sin 8 dO 

The thermal energy equation can be re- written as (see SZ92) 
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where g — d<I>/dr and N is the Brunt- Vaisala (buoyancy) 
frequency. 

The equation for the conservation of magnetic flux can be 
integrated once, and axial symmetry implies that 
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The (^—component of the equation of conservation of mag- 
netic flux: 
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Finally, the solenoidal condition is 
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APPENDIX C: NUMERICAL 
IMPLEMENTATION OF THE EQUATIONS 

Numerical implementation of the equations is done by defin- 
ing the non-dimensional independent variables x — r / Rq 
and p = cos#, and by working with the non-dimensional 
dependent variables 
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The perturbed thermodynamical quantities are also normal- 
ized with 

p=-?- ,f=^ ,p=^- , (C2) 

PO -LO PO 

where po = lg/cm 3 , To = IK, and po = poRqQ 2 . 

The symmetries of the system in the latitudial coordi- 
nate 8 suggest the expansion of the dependent variables onto 
Fourier modes, which are equivalent to Chebishev polynomi- 
als of the variable p since the n— order Chebishev polynomial 
T n is defined as 



T„(p) = cos(n6 



(C3) 



We assume symmetry across the equator for the radial ve- 
locity, the angular velocity, the thermodynamical quantities 
and the latitudinal component of the magnetic field. Anti- 
symmetry across the equator is then required for consistency 
for the latitudinal velocity as well as the radial and toroidal 
components of the magnetic field. In addition, we require 
that the total mass flux across a spherical surface be null. 
This leads to the suggested expansion: 

N d 
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JV 

f(x,p) = ^e„wr 2 „- 2 ( M ), 



B(x,p) = ^B B (x)T2»-j(/i), 

JV 

b(x,p) = (1 - p 2 ) b n (x)T 2n ^ 2 (p) , 

n = l 
JV 

S(x,p) = (1 - p 2 ) ^2 S„(x)Tan-i(n) , 

n=l 
JV 

j(x,p) = (l-p 2 )Yl Mx)Tm-iO*) . (C4) 



where each sum is truncated to retain the first N modes only. 
The equations are expanded using the ansatz l|C4|> . then 
projected back onto the first N Chebishev modes (of rele- 
vant parity) . The quadratic terms are simplified analytically 
using the standard formulae 2T n (p)T m (p) = (T n+m (p) + 
Tn-m(p)) and T„(p) = T- n (p), as well as a variety of oth- 
ers that have been summarized by Garaud (2001, pp. 200- 
204). This procedure yields a combination of 12iV first-order 
ODEs and 3jV algebraic equations for a total of 15jV depen- 
dent variables, namely {ipn}, {v n }, {L n }, {T„} and {5*„} 



© 2007 RAS, MNRAS 000,rTH2T1 



20 P. Garaud & J.-D. Garaud 



(and their respective radial first derivatives), as well as {p n }, 
{p n }, {B n }, {b„} and {J„} for n = 1..N. 

These equations are solved using a versatile variant of 
the standard Newton-Raphson-Kantorovich (NRK) relax- 
ation solver (see for example Press et al. 1 996, chapter 17.3) 
developped by Garaud (2001, pp. 108-111). The added fea- 
ture compared with the standard algorithm permits the so- 
lution of any set of equations that can be written in the 
form 



E 

3=1 



fi(y,x) for i = l.JVv 



(C5) 



where x is the independent variable, y = {y;} is the vector 
containing the iVv dependent variables, and Mij(y,x) and 
/i(y, x) can be any nonlinear function of x and y. One of the 
many advantages of this form is that algebraic equations are 
straightforwardly treated as any other equation by setting 
Mij = for all j and for i corresponding to the relevant 
equation(s). 

As in the standard relaxation algorithm (see Press et 
al. 1996) most of the computational cost (both in terms 
of time and memory) arises from reading and inverting a 
block-tridiagonal matrix containing N x + 1 rows of 3 blockfQ 
(where N x is equal to the number of meshpoints used, typi- 
cally between 2000 and 3000) , and where each block has siz^fl 
N v x N v where N v is the total number of dependent vari- 
ables (N v — 15N, where N is 60-80 for a typical run). Cal- 
culation of the matrix components, followed by its inversion 
using serial partial pivoting can require several hours per it- 
eration on a high-end desktop. In addition, memory becomes 
an issue since the standard algorithm typically requires the 
storage of N x double-precision arrays of size N v X JV v /2 be- 
fore back-substitution can proceed. For this reason, we have 
developped a parallel version of the NRK solver. 



APPENDIX D: PARALLEL NRK ALGORITHM 

A serial, simple and efficient way of inverting the typi- 
cal block-tridiagonal linear systems arising from two-point 
boundary value problems is described by Press et al. (1996). 
Given the matrix structure 
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the algorithm consists in diagonalizing the first top M (mid- 
dle) block using Gauss- Jordan elimination with partial im- 
plicit pivoting, and storing the resulting modified R block 
and right-hand side for later back-substitution. Moving onto 
the second block-row, we zero the L block, diagonalize the 



7 except for the boundary blocks, of course. 

8 again, except for the boundary blocks, which are smaller. 



M block and store the modified R block (and modified right- 
hand-side vector) , and so forth. The resulting matrix struc- 
ture is 
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(D2) 



where denotes a "zeroed" block, / is the identity, and S is 
a stored block. The last step involves sweeping the matrix 
back from bottom to top to back-substitute for the unknown 
variables. This method requires minimal storage: firstly, each 
block-row only needs to be read just before being processed 
and secondly, the backsubstitution step only requires knowl- 
edge of the blocks S and modified right-hand-side vector. 

Inverting the same linear system in parallel is not obvi- 
ous since the previous algorithm is inherently serial. A stan- 
dard parallel method used for block tridiagonal matrices is 
cyclic reduction (see Golub & Ortega, 1993). In this case, 
the number of variables is successively halved by repeating 
the following reduction algorithm: diagonalize all M-blocks 
in odd block-rows by partial implicit pivoting, then zero the 
R and L matrices directly above and below by Gaussian 
elimination. The resulting matrix after the first reduction 
is: 
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(D3) 



where S denotes a block that has been modified and needs to 
be stored for later back-substitution, and X denotes blocks 
that have been modified and will be used in the next re- 
duction step. Note that all of the remaining X blocks (in 
the odd block-rows) form a new block-tridiagonal system 
for roughly half the number of variables (more precisely, 
if the initial number of block-rows is 2 k — 1 then the re- 
duction step reduces it to 2 fc ~ 1 — 1). The algorithm can be 
repeated until only 1 block is left and finally inverted, at 
which point backsubstitution can begin. The main advan- 
tage of the method is its inherently parallelizable nature, 
since each process can be instructed to reduce a number of 
blocks more-or-less independently of the others; communica- 
tion between processes is minimal until the final steps where 
the number of remaining block-rows is equal to the number 
of processes. Unfortunately, the stability of this cyclic reduc- 
tion algorithm is much weaker that the stability of the serial 
algorithm (since the matrix cannot be globally pivoted), and 
was found to be unreliable for our purpose. 

An alternative parallel algorithm was constructed, 
loosely based on the work of Wright (1992). The block-rows 
are first equally distributed between the processes. A first 
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sweeping reduction akin to the serial algorithm is performed 
starting from the top block-row down to the second-to-last 
block-row in each proces^f] to transform the matrix to 
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(assuming in the case shown here that 12 block-rows are 
equally distributed between 3 processes). Again, by conven- 
tion all S blocks denote blocks that need to be stored for 
backsubstitution, while X blocks are blocks that will be fur- 
ther modified. Until this point, work in each process can be 
done without communication with other processes. 

At the end of the sweeping step, each process sends its 
very last block-row to the following process. Thus the blocks 
are redistributed as 
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Each process can then continue eliminating undesirable vari- 
ables from its top block-row by Gaussian elimination with 
the successive block-rows below, until the following form is 
achieved: 
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Figure Dl. Scaling properties of the parallel NRK algorithm 
RELAX for the simulations presented in this paper. The stars 
correspond to a simulation with 60 modes (i.e. 900 equations), 
and 1000 mcshpoints. The diamonds correspond to the same sim- 
ulation but with 2000 mcshpoints. 



At this point, the top-rows in each process can be com- 
bined into a much-reduced block tri-diagonal system. So- 
lution from here on proceeds with a cyclic reduction of 
the remaining blocks across processes, followed by a back- 
substitution step. 

The operation count of this algorithm is the same as 
that of the standard cyclic reduction. The advantages of 
this algorithm over the standard cyclic reduction are two- 
fold. Firstly, it is found to be much more stable (Wright, 
1992). Secondly, it is far more versatile, since the standard 
algorithm only performs ideally with a number of block-rows 
equal to 2 fe + 1 (for any integer k) , which seriously constrains 
the number of meshpoints one is allowed to select. On the 
other hand, in this algorithm it is always possible to balance 
the load (for any number of meshpoints) in such a way that 
any process contains at most one more block-row than the 
average. 

The scalability of the algorithm is naturally excellent 
until the number of processes approaches a tenth of the num- 
ber of meshpoints. This is illustrated in Figure iDll 

Finally, we note that the user-interface of the com- 
plete associated parallel version of the Newton-Raphson- 
Kantorovich algorithm (RELAX) is entirely identical to that 
of the widely used serial NRK algorithm written by Gough 
& Moore. Transfer between the parallel and serial version 
should be transparent to the user. The RELAX algorithm 
is freely available upon request from P. Garaud. 

This paper has been typeset from a H3X/ LTffjX file prepared 
by the author. 



except in the last process, which is swept all the way 
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